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In this paper, an analytical solution of alloy solidification problem is presented. 
We develop a special method to obtain an exact analytical solution for mushy zone 
problem. The main key of this method is a requirement that thermal diffusivity in the 
mushy zone to be constant. From such condition we obtain an ordinary differential 
equation for liquid fraction function. Thus the method can be examine as "a model" 
to achive analytical solution of some unrealistic problems. 

An example of solutions is presented: the noneutectic titanium-based alloy 
solidification. We provide the comparison of numerical simulation results with obtained 
exact solutions. It shown that very simple apparent capacity-based numerical scheme 
is provided a good agreement with exact positions of the solidus and liquidus isoterms, 
and with temperature profiles also. 

Finally, some extensions of the method are outlined. 



1. Introduction. 

A general methodolgy of achieving analytical solutions of the alloy solidification problem 
is presented in this manuscript. There is an analytical solution for pure substance 
(Stefan's problem) and few analytical and semi-analytical solutions for alloys [1, 2, 3]. 
We suggest a general methodology which can provide wide range solutions to test 
different numerical schemes [1, 2]. 

We consider the case when physical properties $ (density, heat capacity, or heat 
conductivity) in solid and in liquid are constant. Within mushy zone these properties 
and the enthalpy depend on temperature as follows (i.e. obey the lever rule): 



where $s and $/ are properties in solid and liquid, respectively, A(T) is volumetric liquid 
fraction. Then we rewrite heat transfer equation in the full enthalpy term H{T). The 
key idea of the present work is the mushy heat diffusivity requirement to be constant. 



Abstract. 



$(r) = [i-A(T)]$, + A(r)<i>z, 



(1) 



a(T) 



Gsi = const, 



(2) 



dH{T) 
dT 
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where k{T) is heat conductivity. From this condition we can find hquid fraction A(T) by 
means of which we are able to hnearize an initial energy conservation equation. Thus, 
the following methodology is 

(i) To rewrite of the heat equation in the full enthalpy term. 

(ii) To require of the thermal diffusivity to be constant in the solid, mushy and liquid 
zone. 

(iii) Condition a{T) = agi — const is ordinary differential equation for liquid fraction 
A — A(T). Additionaly we require A(T;) = 1. 

(iv) To solve this ODE and find A = A(r, asi). 

(v) To impose additional condition \(Ts,asi) = Xq. For Aq = we have noneutectic 
alloy, and for Aq 7^ - eutectic. From this condition we find agi- 

(vi) Now we have the heat equation with constant-peace coefficients and we can it solve 
easy. 

It needs to note, that this problem cannot be solved with well defined (predefined) A(T) 
function, instead of the function A = A(T) is determined from linearisation conditions. 



2. The linearisation of the heat equation. 

We will solve an energy conservation equation 



dH 

— = div {k{T)grad T) , (3) 
where full enthalpy H is 

T T 

H{T) = p{T) [1 - A(T)] j CsiOdC + P{T)X{T) j Ci{C)dC + p{T)X{T)L, (4) 



where Cg, Ci is specific heat in solid and in liquid, L is latent heat of fusion, p = Ps = Pi 
is density, which all are constants. We express the heat conductivity in the "mixture" 
form 

k{T)^[1-X{T)]ks + X{T)ku (5) 

where Kg and ki are constant heat conductivity in solid and liquid, respectively. 
Taking into account the expression 

gradH 

gradT^ (6) 

dT 

the Eq. (3) can be rewritten in the general form 

— = div{a{T)gradH), (7) 
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where a{H) is thermal diffusivity, which defined as 

dT 

If K,{T) and dH{T)/dT depend on temperature arbitrary manner then Equation (7) is 
nonhnear. To achieve an anahtical solution we need to require the thermal diffusivity 
to be constant in all regions (solid, mushy and liquid). 





OLs — const 


for T <Ts 




a(T) = < 


asi — const 


for Ts<T <Ti 


(9) 




ai — const 


for T>Ti 





In our case = KgpCs and ai = Ki/pCi are constant by definition. For the derivation 
of mushy enthalpy (the apparent capacity x density) we get: 

= [1 - A] + XQ + [(Q - C^)T + L] (10) 

then from the mushy part of Eq. (9) we obtain an ordinal differential equation for 

A(r) 

[l+pT]^ + aX{T) + b^O, (11) 
where we denote 

_ asipjCi - Cs) - {ki - Ks) . s 

a — , \i.Z) 

ttsipL 



O-slpCs — Us 

ttsipL 



(13) 



Ci — Co 

(14) 



We require 



KTi) = 1, (15) 
where T; is a liquidus temperature. Solution of Eqs. (11) and (15) is 

A(r, = -^ + ^('i±P^V. (16) 

It needs to determine an additional condition for X{T) function, namely to define the 
liquid fraction value at solidus temperature 

{0 for noneutectic alloy 
Ao for eutectic alloy 

To obtain the analytical solution of Eq. (7) we need to solve Eq.(17) to find root a^/. 
Then we need to solve Eq. (7) with suitable initial and boundary conditions. 
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CsT for T<Ts 

CsT + \{T)[{Ci-Cs)T + L] iorT,<T<Ti . 
QT + L iorT>Ti 



H{T) 



In future we need the value dT/dH: 
1 



dT_ 
'dH 



1 



dX 
dT 



for T<Ts 
for Ts<T <Ti 

for T>Ti 



(19) 



[l-X]Cs + XCi + [{Ci-CsT + L)] 
1 

It should be note that some expressions with dX/dT may be written in more 
simplified form versus A, for example: 

dX aX + b aX + b 

dT^~ 



1+pT (Ci-Cs)T + L 



(20) 



and the combination 



m-C,)T + L]^^-{aX + b)L. 



3. An analytical solution for enthalpy. 

We will examine the simple problem 



(21) 



dH d ( ,„.dH\ 

H{t = Q) = H{Tir^u), 



H\ 



z=0 



Hout 

--H, 



init 



H(Tout), 
= H(Tinit) 



(22) 

(23) 
(24) 
(25) 



The solution of these equations with constant-piece function a{H) can be easy find [4]. 
To solve this equation we divide whole region [0, oo) into three subintervals [0,Xs), 
[Xs,X;] and {Xi.oo) {Xg and Xi are solidus and liquidus positions, respectively). 
Moreover, we assume that (the similarity solution): 



Xs{t)^ks\ft, Xi{t)^kiVt, 



(26) 



where kg and ki are constants. Solutions on the subintervals are: 
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H{x,t) 



{Hi H,)n 



2v <>.«/' 



erf 



H{x,t) = Hinit - {Hinit - Hi 



where we defined 



erfc 



f^, xe{Xi,oo), 



, xe[X,,Xi], (28) 



(29) 



H, = H{Ts) = pCsTs, 
Hi^H{Ti) = p{CiTi + L). 



(30) 
(31) 



By using the two conditions at the interfaces (the first one from which is Stefan's 
condition at the sohdus (eutectic) point): 



dH 




dH 


ox 


x=Xs-0 


ox 


dH 




dH 


ox 


x=Xi-0 


ox 



+ pXoL 



x=Xs+0 



x=Xi+0 



dt ' 



(32) 
(33) 



we derive the following two equations from which to evaluate ks and ki: 
^s{Hs - H^t)exp (-^) ^i{Hi - Hs)exp [S]) _ 



eri 



/(*) 

oTiiHi - Hs)exp ^ 4^ 



Oii{Hinit - Hi)exp 



er 



(34) 



(35) 



4. What we can get from the exact solution? 

Usualy wc have numerical scheme which gives us the temperature, but not enthalpy. 
Below wc write down formulas for temperature evaluation versus enthalpy and some 
other parameters. 



4-1- Solidus and liquidus velocities. 

From Eq. (26) we get the front velocities 

dXi{t) 



vi(t) = 



dt 



ki 
2x/t' 



Vs{t) = 



dXs{t) ks 



dt 



2y/t' 



(36) 
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Prom Eq. (18) and Eqs. (27) - (29) we can easy to find: 
In the solid {x < kgVt): 



T{x,t) 



pCs 



Hout + {Hg — Hf 



erf 



2vW 



out I 



er 



(37) 



In the mushy zone {ksVt < x < kiy/t) we need to solve nonlinear equation to get 
T^T(x,t): 



C,T + X{T) [{Ci-Cs)T + L] 

_l {Hi-H,)erf, ,^, 
P 



{2^) - Hierf (^) 



(^fc) - ^"-f (27^ 



T{x,t)^ 



In the liquid {x > ki^/t): 
1 



pCi 



Hinit ~ {Hinit ~ Hi 



L 



4.3. Temperature gradients. 
From equation 



dT dT dH 



dx dH dx 

end from Eq. (19) easy to achive the expressions for temperature gradients: 

r 1 OH 



dT 
dx 



pCs dx 
1 



1 



dH 



p [1 - A] + XCi + [(Q - Cs)T + L] § ax 
1 dH 
. pCi dx 

For enthalpy gradients from Eqs. (27) - (29) we get 



r<r. 



Ts<T<Ti 



T>Ti 



dH 
dx 



Hs — H, 



out 



x2 



er 



Hi — Hs 



Hinit ~ Hi 



e "'"i' 



- er/c (^) ■ x/^ 



H<H. 



Hs<H< Hi 



H>Hi 



(38) 



(39) 



(40) 



(41) 
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dT 
dx 



1 H. - H. 



out 



(Hi - Hs)/ (erf (^) - erf 



[l-X)Cs + \Ci + [{Ci - Cs)T + L]§ 



X < ksVi 

ksVi <x< kiy/i (42) 



init 



H, 



pCi erfc 



x2 
\/T^OLlt 



X 



> kVi 



There is a very interesting parameters as temperature gradient in liquid phase at 
the liquidus point Gi. For it we can write down 



Giit) 



dT{x, t) 



dt 



Hinit ~ Hi 



=Xi+o P^i erfc 



(43) 



This parameter controls the type of solidification microstructure. 



4- 4- Cooling rate. 

As we can see Giit) ~ 1/ \/t. Liquidus velocity (36) varies with time alse as vi(t) ~ 1/ ^/t, 



then cooling rate given by GiVi ^ 1/t. It is easy to show that is so. The cooling rate is 
defined as 

dT{x, t) dT dH 



T{x,t) 



dt 



dT 
'di 



— < 



Hs — H, 



out 



dH dt 

asXe *°:st 



(Hi - Hs)/ {erf (^) - erf (^)) 



asixe ^"si^ 



(1 - A)C, + XCi + [{Ci - Cs)T + L] §] 20F(a,ii)3/2 



init 



H, 



aixe 



^2 
4a; t 



pCi erfc{^) 2yi(a,t)3/2 
Thus, cooling rate at the liquidus point is given by 



(44) 
X < ksVi 

ks\/i < X < ki\/t 

X > kiy/i 



T,^T 



H., 



init 



Hi kic ""i 



x=Xi 



+0 pCi erfc[^) 2^/Wait t' 



(46) 



The value Ti is very important, because it defines secondary arm dendrite spacing 



[7]. The another important expression is ^^^W; which defines primary arm dendrite 
spacing [8]. From Eqs. (36) and (43) we can show that primary arm spacing varies versus 
time like ~ t^^^. However, as it's very known, after some critical gradient and velocity 
at the liquidus point will be take a place columnar to eqiaxed transition [9]. 
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4.5. Local solidification time 

For directional solidificcation, the local solidification time tis can be estimated from the 
following equation: 



tls{x) 



1 1 



(47) 



where quadratic increasing with x oitis we have, because the mushy zone lehgth increases 
versus x and sohdus/liquidus velocities decrease. The local sohdification time controls 
the some segragation processes in the mushy zone. 

5. Numerical scheme. 

We will solve heat transfer equation 



(48) 



which concerns with Eq. (23). Here k{T) edfined by Eq. (5) and mushy heat capacity 
(so-called apparent capacity): 

d\{T) 



C{T) = [1 - A(T)] Cs + \{T)Ci + [{Ci - Cs)T + L] 



dT 



(49) 



The first, we draw the grid with spatial step h = Xj+i — Xi, where i = 0,1, ...,N. 
The second, integrating the Eq. (48) over the x e [(i — ^)h, {i + ^)h] we can write heat 
balance equation as follows 



Ti+i — Ti 



h 



x={i-^)h 

The left part of this equation we express as 



' 2 



Ti — Tj-i 
1 



(50) 



Ot T 

where n is time index, r is time step. Thus the Eq. (48) we can write down in the 

discrete form as 



' 2 



where [10] 



(51) 



^j-i = 

2 + Ki 



2 f^/j K/'. 



+1 



2 Ki + Ki+i 

We note, that Cj etc are calculated at time tn, i.e. 



Q = C{T^), = k{TI^ etc. 
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Boundary conditions are 

To = Tout, Tn = Tinit- (52) 

Because the Eq. (51) is tri-diagonal linear system, then its solving is trivial and we do 
not discuss this issue here. 

6. Binary alloy solidification: numerical treatment. 

In this section we consider the solidification of noneutectic titanium-based alloy, which 
we can treat as pseudo-binary alloy. Physical properties of titanium alloy VT3-1 (Ti- 
6.5Al-2.5Mo-1.5Cr-0.5Fe-0.3Si) are present in the Table 1. These parameters we are used 
for numerical simulation of liquid pool profiles during vacuum arc remelting process [5] . 



Table 1. Properties 
of the VT3-1 alloy. 

Parameter Value 



Cs 600 J/kgK 

Ci 1200 J/kg K 

Ks 10 W/m K 

Ki 35 W/m K 

p 4500 kg/m^ 

L 3.55 X 10^ J/kg 

1550 °C 

Ti 1620 °C 

Tm 1668 °C 



A solution of the Eq. (17) with Aq = is agi = 2.26891 x 10 m^/s. Figure 1 shows the 
temperature dependence of the liquid fraction. Additionaly Figure 1 shows the function 

XtiT) = 1 - • (53) 

which we used for VT3-1 alloy [5]. The difference between Ay and Xt{T) is small, then 
we have nearly realistic problem, li g(T) approximated with a power function [6] 



A^ 



^Ti-T,, 

then we get n ^ 1.5. To test very simple numerical apparent capacity-based method 
we carried out simulations with following parameters: = 3.7037 x 10^® m^/s, 
ai = 6.48148 x 10"^ m'^/s, Tout = 800 °C, Tinu = 1650 °C, Hout = 2.16 x 10^ J/m^, 
Hinit = 10.5057 X 10^ J/m^, = 4.185 x 10^ J/m^, Hi = 10.3455 x 10^ J/m^. 

Solutions of Eqs. (34)-(35) are ks = 0.00134109 m/s^/^ k = 0.00206009 m/s^s. 

The numerical model parameters are chosen as: length of domain d — 0.5 m, 
nodes number N — 500, time step r = 0.1 s. A numerical model can provide excellent 
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Figure 1. Liquid fraction versus temperature for noneutectic alloy. On the Figure: 
analytical - for X{T), VT3-1 for At(T). 



agreement with obtained analytical solution. The results obtained are in a Figure 2: 
movement of both the solidus and the liquidus front (26). We used linear interpolation 
between Tj and Tj+i for estimating position of the fronts, whereas Tg^i G [Tj,Tj_|_i]. 




Figure 2. Comparison of the apparent capacity-based numerical method with the 
analytical solution as applied to the solidification of VT3-1 titanium alloy. Figure: 
solidus and liquidus position versus time. 



The errors in the positions of solidus/liquidus 

£x[t) = 7-T iUU/o 

-^exact ) 



are presented in the Figure 5. Moreover Figure 4 shows the temperature profiles after 
20 and 500 seconds under the same numerical conditions. 
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Figure 3. Solidus/liquidus front position errors versus time. 
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Figure 4. Temperature curves at 20 and 500 seconds. 
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Figure 5. Temperature curves errors at 20 and 500 seconds. 



Towards analytical solutions of the alloy solidification problem. 



12 



Figure 5 shows temperaure profiles errors defined by 

St{x) = — 100%. 

exact y-'^ ) 

We would like to underline that the purpose of this work is to achieve the exact 
analytical solution on alloy solidification. Due to this, advantages and disatvantages 
different numerical algorithms can be done in future. Due to this, wc don't study 
the process of solidification of Ti-6.5Al-2.5Mo-l.5Cr-0.5Fe-0.3Si alloy, but only use the 
thermo-physical properties of this alloy to show as the model works. 



7. Conclusions. 



In this paper, analytical solution of alloy solidification problem is presented. We 
developed a special method to obtain an exact analytical solution for mushy zone 
problem. The main requirement of the method is thermal diffusivity to be constant 
in the mushy zone. Due to such condition ordinary differential equation for liquid 
fraction function is achieved. Thus the present method can be examined as "a model" 
way to get analytical solution of some unrealistic problems. 

An example of solutions is given - the noneutectic titanium-based alloy 
solidification. We provide the comparison of the simple numerical simulation results 
with obtained exact solutions. We show that very our numerical apparent capacity-based 
scheme provides a good agreement with exact solutions for solidus/liquidus position and 
for temperatures profiles in different moments of solidification time. 

Once again we would like to underline that the main goal of this paper to provide 
the benchmark for binary alloy solidification problem, but not in the analysis of used 
numerical scheme. 

If predefined A(T) function is to be examined, we can use another suggestions. For 
example, we can require to heat conductivity (from experiment, e.g.) to be proportional 
to the apparent capacity, i.e. 

k{T) = asip [Cs + {Ci - Cs)X{T) + [{Ci - Cs)T + L] j . 

Or, for the second example, we require to apparent capacity (from experiment) to be 
proportional to mushy heat conductivity, i.e. 

dH{T) ^ Ks + {ki - Ks)X{T) 
dT Gsi 

Moreover, we may use the Backlund's transformation [11] to make mushy heat equation 
linearisation. In this case we get nonlinear condition 

H\T)X{T) 



dH{T) 

dT 



const. 



These linearization methods will provide us with some additional analytical solutions of 
alloy solidification problem. 
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